#Clear R
rm(list=ls())
zero1 <- function(x, minx=NA, maxx=NA){
  res <- NA
  if(is.na(minx)) res <- (x - min(x,na.rm=T))/(max(x,na.rm=T) -min(x,na.rm=T))
  if(!is.na(minx)) res <- (x - minx)/(maxx -minx)
  res
}

load("Kam_data.RData")

#function to install packages if they don't exist
ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}

# usage
packages <- c("ggplot2", "psych","interplot","stargazer","dplyr","stringr","tidyr", "lavaan")
ipak(packages)

###############
#Descriptive statistics in paper
###############
mean(data$DV_irradiation,na.rm=T)
sd(data$DV_irradiation, na.rm=T)
alpha_dv<-psych::alpha(data.frame(data$SupportIrradiation_rec, data$CostsFoodIrradiation_rec, data$IrradiationGood_rec))
summary(alpha_dv)


############################################
#
# Figure 1
#
#############################################

#model ANES
summary(NfCANES<-(lm(data$DV_irradiation~InParty*NFC_ANES +NoCues*NFC_ANES,data=data)))
#make a grpah
df1<-(interplot(m = NfCANES, var1 = "InParty", var2 = "NFC_ANES"))
df1<-(data.frame(df1$data))
df1$index<-"ANES 2-item"

#modle Bullock
summary(NfCBullock<-(lm(data$DV_irradiation~InParty*NFC_Bullock +NoCues*NFC_Bullock,data=data)))
#make a graph
df2<-(interplot(m = NfCBullock, var1 = "InParty", var2 = "NFC_Bullock"))
df2<-(data.frame(df2$data))
df2$index<-"Bullock 6-item"

#model full 18-item
summary(NFC18<-(lm(data$DV_irradiation~InParty*NFC18 +NoCues*NFC18,data=data)))
#make a graph
df3<-(interplot(m = NFC18, var1 = "InParty", var2 = "NFC18"))
df3<-(data.frame(df3$data))
df3$index<-"Cacioppo et al. 18-item"

#create plot
toggplot<-(rbind(df1,df2,df3))
ggplot(toggplot,aes(x=fake,y=coef1))+geom_line()+geom_ribbon(aes(ymin=ub, ymax=lb), alpha=0.5)+facet_wrap(~index)+ geom_hline(yintercept=0, linetype = "dashed") + ylab('Marginal effect of In-Party Cue on Support for Policy') +xlab('Need for Cognition')+theme_bw()
ggsave("~/Dropbox/Apps/ShareLaTeX/TIPI Validation Study/Figures/Kam_Figure3_final.pdf",width=10,height=8)

#Extract results for models
stargazer(NfCANES, NfCBullock, NFC18, title="Moderation of Party Cues", align=TRUE, 
          dep.var.labels=c("Support for Ban of Food Irradiation"), 
          covariate.labels=c("In-Party Cue", "Need for Cognition", "No-Party Cue", "In-Paty Cue X Need for Congition", "No Cue X Need for Cognition"), 
          omit.stat=c("LL","ser","f", "adj.rsq"), 
          star.cutoffs=c(0.05), 
          notes = "OLS Regession models; *p<0.05", 
          notes.append = FALSE, 
          no.space=TRUE)

##############################
#
#Figure 2
#
##############################
nc <- data
## Put NC items into own dataframe
ncitems <- dplyr::select(nc,c("NfC_1","NfC_2","NfC_3_rec","NfC_4_rec","NfC_5_rec","NfC_6","NfC_7_rec","NfC_8_rec","NfC_9_rec","NfC_10","NfC_11","NfC_12_rec","NfC_13","NfC_14","NfC_15","NfC_16_rec","NfC_17_rec","NfC_18"))
names(ncitems)

CF18<-
'NfC_full18 =~ NfC_1 + NfC_2+NfC_3_rec+NfC_4_rec+NfC_5_rec+NfC_6+NfC_7_rec+NfC_8_rec+NfC_9_rec+NfC_10+NfC_11+NfC_12_rec+NfC_13+NfC_14+NfC_15+NfC_16_rec+NfC_17_rec+NfC_18
# fix variance of speed factor
NfC_full18 ~~ 1*NfC_full18'

#assess fit
fit<-cfa(CF18,data=ncitems)


## Create Empty List to stick output into  k = the number of combinations when we have between 2 and 17 items in a set
k <- 0
for(i in 2:17){
k <- k + ncol(combn(1:18, i))
}

nclist <- list()
# empty matrix to stick alpha scores
alphacoefs <- data.frame(matrix(nrow=k,ncol=4))
h=1
for(j in 2:17){
  # create possible combinations of item numbers of length j
  combos <- combn(1:18, j)
  # empty matrix to stick individual means made of of items from combinations set in combos
  
  two <- matrix(nrow=nrow(nc),ncol=ncol(combos))
  for(i in 1:ncol(combos)){
    two[,i] <-   rowMeans(ncitems[,combos[,i]])
    alphacoefs[h,] <-   psych::alpha(ncitems[,combos[,i]],check.keys = TRUE)$total[1:4]
    h=  h+1
  }
  ## m
  
  lmcoefs <- matrix(nrow=ncol(combos),ncol=4)
  for(i in 1:ncol(two)){
    lmcoefs[i,]   <-  summary(lm(zero1(nc$DV_irradiation)~zero1(two[,i])*nc$InParty+zero1(two[,i])*nc$NoCues))$coefficients[5,]
  }
  print(j)
  nclist[[j]] <- lmcoefs
}

b <- 0
items <- list()
for(j in 2:17){
  b=b+1 
  combos <- combn(1:18, j)
  for(i in 1:ncol(combos)){
    items[[b]] <- combos
  }}

itemrow <- lapply(items,nrow)
itemcol <- lapply(items,ncol)

collapsepastecolumns <- function(x)apply(x,2,paste,collapse=",")
rr <- lapply(collapsepastecolumns,X = items)

numbero <-   unlist(itemcol)
emptylist <- list()  
for(i in 2:17){
  emptylist[[i-1]] <- rep(i,numbero[[i-1]])
}



coefficients <- do.call("rbind",nclist)
itemstogether <- unlist(rr)
numall <- unlist(emptylist)


library(ggplot2)
both <- data.frame(coefficients,alphacoefs,numall,itemstogether)
head(alphacoefs)
eighteen <- rowMeans(ncitems,na.rm=T)
bottomrow <- summary(lm(zero1(nc$DV_irradiation)~zero1(eighteen)*nc$InParty+zero1(eighteen)*nc$NoCues))$coefficients[5,]
names(bottomrow) <- names(both)[1:4]

oo <- rbind(both,NA)
oo[nrow(oo),1:4] <- bottomrow
oo[nrow(oo),]$numall <- 18


for(i in 1:18){
  bottomrow <- summary(lm(zero1(nc$DV_irradiation)~zero1(ncitems[,i])*nc$InParty+zero1(eighteen)*nc$NoCues))$coefficients[5,]
  oo <- rbind(oo,NA)
  oo[nrow(oo),1:4] <- bottomrow
  oo[nrow(oo),]$numall <- 1
}

ggplot(oo,aes(y=X1.1,x=as.factor(numall)))+geom_boxplot()+theme_bw()+xlab("Number of items")+ylab("Coefficient")

ggplot(subset(oo,numall!=1 & numall!=18),aes(y=X1.1,x=X1))+geom_point()+theme_bw()+xlab("Cronbach's Alpha")+geom_smooth(method="lm")+ylab("b Coefficient")+facet_wrap(~numall)

## split up which items went into which coefficient
na.omit(oo) %>% separate(itemstogether, into = paste("V", 1:18, sep = "_")) -> both
itemscores <- inspect(fit,what="std")$lambda
itemsforfactor <- dplyr::select(both,starts_with("V_"))
for(i in 1:18){
  itemsforfactor[itemsforfactor==i] <- itemscores[i]
}

itemsforfactor <- mapply(as.numeric,itemsforfactor)
factormeans <- rowMeans(itemsforfactor,na.rm=T)
both$numall_1 <- factor(paste(both$numall,"Items"))
both$numall_1 <-  factor(both$numall_1,levels=levels(both$numall_1)[c(9:16,1:8)])
all <- data.frame(both,factormeans=factormeans)

ggplot(all,aes(y=X1.1,x=factormeans))+geom_point()+theme_bw()+xlab("Mean Factor Score of Items")+geom_smooth(method="lm")+ylab("b Coefficient")+facet_wrap(~numall_1)

